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Abstract 



We present a cliaracteristic algorithm for computing the perturbation of 
a Schwarzschild spacetime by means of solving the Teukolsky equation. We 
implement the algorithm as a characteristic evolution code and apply it to 
compute the advanced solution to a black hole collision in the close approxi- 
mation. The code successfully tracks the initial burst and quasinormal decay 
of a black hole perturbation through 10 orders of magnitude and tracks the 
final power law decay through an additional 6 orders of magnitude. Deter- 
mination of the advanced solution, in which ingoing radiation is absorbed by 
the black hole but no outgoing radiation is emitted, is the first stage of a two 
stage approach to determining the retarded solution, which provides the close 
approximation waveform with the physically appropriate boundary condition 
of no ingoing radiation. 

PACS number (s): 04.20.Ex, 04.25.Dm, 04.25.Nx, 04.70.Bw 
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Typeset using REVT^ 



I. INTRODUCTION 



In this work, we present the advanced solution for a perturbation of a Schwarzschild 
background describing the head-on coUision of black holes in the close approximation where 
the merger takes place in the far past. We compute the solution by means of a characteristic 
evolution of the Weyl tensor, as governed by the Teukolsky equation [|l],0. The advanced 
solution corresponds to Stage I of a new two stage approach to the vacuum binary black 
hole problem 0. In subsequent work (Stage II), we will use the results of this first stage to 
compute the physically more relevant retarded solution. This perturbative solution will in 
turn provide a valuable reference point for the physical understanding of a fully nonlinear 
treatment of binary black holes being pursued by a similar two stage strategy IQ, a compu- 
tationally feasible problem using an existing characteristic code |^. The perturbative results 
also provide a new perspective on the physical picture previously obtained by applying the 
Cauchy problem to the close approximation especially with regard to global issues which 
have not yet been explored in the Cauchy formulation. 

In the characteristic formulation of this problem, the advanced solution is simpler to 
compute than the retarded solution because of the global relationship between the null 
hypersurfaces on which boundary information is known. One of these hypersurfaces is the 
black hole event horizon Ti"*" whose perturbation must correspond to the close approximation 
of a binary black hole. In the retarded problem, the other null hypersurface (in a conformally 
compactified description) is past null infinity where the incoming radiation must vanish. 
Because and T~ are disjoint, there is no direct way to use data on those two hypersurfaces 
to pose a characteristic initial value problem. 

However, in the advanced problem, it is at future null infinity that the radiation is 
required to vanish. Since Ti^ and formally intersect at future time infinity /+, they 
can be used to pose a characteristic initial value problem to evolve backward in time and 
compute the exterior region of spacetime. Potential difficulties in dealing with I'^ are avoided 
by posing the no outgoing radiation condition on an ingoing null hypersurface which 
approximates X"*" and intersects T^"*" at a late time when the perturbation of the black hole 
has effectively died out. These data on Ti^ and J"*" then constitute a standard double-null 
initial value problem ||7|-[10| for the exterior spacetime, in which ingoing radiation is absorbed 
by the black hole but there is no outgoing radiation. 

This advanced solution provides the radiation incident from X^. In Stage II of the 
approach, this ingoing radiation will be used to generate a "source free" advanced minus 
retarded solution. A purely retarded solution can then be produced by superposition with 
the Stage I solution. Although we do not address Stage II in this paper, we will discuss 
the role of time reflection symmetry in the perturbation equations, which simplifies the 
technical details in carrying out the superposition. From a time reversed point of view, the 
Stage I solution is equivalent to the retarded solution for a "head-on" white hole fission, 
with the physically relevant boundary conditions that radiation is emitted but that there 
is no incoming radiation from X~. It is convenient here to formulate the Stage I results in 
terms of such a white hole fission since the characteristic evolution then takes the standard 
form of being carried out forward in retarded time. 

The close approximation has been extremely useful for testing fully nonlinear Cauchy 
evolution codes. The results of numerical Cauchy evolution and close-limit perturbative 
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theory are in excellent agreement in the appropriate regime, giving great confidence in both 
approaches Furthermore, the perturbative approach provides an important tool for 



the interpretation and physical understanding of those results . 

Clearly, this vital synergism between numerical and perturbative approaches should 
also extend to characteristic evolution. However, in all perturbative studies performed to 
date, the background geometry has either been the Schwarzschild spacetime in standard 
Schwarzschild coordinates or the Kerr spacetime in Boyer-Lindquist coordinates. These co- 
ordinate systems are appropriate for comparison with results from nonlinear Cauchy evolu- 
tion but, to our knowledge, there does not yet exist a treatment of the close approximation 
in terms of null coordinates appropriate for the comparison with nonlinear characteristic 
evolution. 

This work provides such a framework. The methods and results presented here are 
expected to have direct bearing on the study of binary black holes presently underway using 
a fully nonlinear characteristic code. Characteristic evolution has been totally successful 
in evolving 3-dimensional single black hole spacetimes for effectively infinite times (t ~ 
60, DOOM in terms of black hole mass) . Although it is not yet known to what extent 
the characteristic approach can handle the inspiral and merger of binary black holes, it is 
clear that the limitations are due to difficulties in treating caustics and not due to high 
nonlinearity. 

The characteristic data which has been obtained for the nonlinear description of a binary 
black hole spacetime generate close approximation data for a perturbative solution. We 
present here a numerical code to evolve such data as a perturbation of a Schwarzschild 
spacetime in null coordinates. Fortunately for our purposes, the perturbative formalism 
due to Teukolsky is amenable to a reasonably straightforward change of background 



coordinates, as observed in Ref. IT^ 



The Teukolsky equation is based upon the decomposition of the Einstein equations and 
Bianchi identities in terms of a conveniently chosen complex null tetrad, as carried out by 



Newman and Penrose in the early 1960 's. The Newman-Penrose formalism allowed 
Teukolsky to construct a single master wave equation for the perturbations of the Kerr 
metric in terms of the Weyl curvature components 1^4^ (describing outgoing radiation) or 
ipQ (describing ingoing radiation). The Teukolsky formalism provides a completely gauge 
invariant spherical harmonic multipole decomposition for both even and odd parity pertur- 
bations in terms of radial wave equations. For a Kerr black hole with angular momentum, 
there is no similar multipole decomposition of metric perturbations in the time domain (as 
opposed to the frequency domain of Fourier modes). In the non- rotating case, the Teukolsky 
equation reduces to the so-called Bardeen-Press equation [T^ 



Since the 1970's the Teukolsky equation for the first order perturbations around a rotat- 
ing black hole has been Fourier transformed and integrated in the frequency domain for a 
variety of situations where initial data played no role |T7| , p!8[] . In order to avoid the impor- 



tant but difficult problem of prescribing physically appropriate initial data for that equation, 
the computation of gravitational radiation has been restricted to the cases of unbounded 
particle trajectories or circular motion. The first evolution code to integrate the Teukolsky 
equation in the time domain, in Boyer-Lindquist coordinates, was recently developed 
and successfully tested . In order to incorporate initial data representing realistic astro- 



physical initial data for the late stage of binary black hole coalescence, 3 + 1 expressions 



3 



connecting ijj^ and its time derivative to Cauchy data (3-metric hij and extrinsic curvature 
Kij) satisfying the Hamiltonian and momentum constraints 

have been worked out exphcitly ||2l|j20| . |2^ . |23 



In Sec. fT|, we specify a null background tetrad suitable to re-express the Teukolsky 
equations for ipQ and ipi in null coordinates appropriate for characteristic evolution. We also 
discuss various global aspects of these equations which are important for numerical evolution. 



In Sec. pl| , we discuss data for the Teukolsky equation. We present null data for linearized 
Robinson- Trautman solutions, which provide an analytic check on numerical accuracy, and 
null data for the close approximation to a white hole fission. In Sec. |V|, we discuss the 
numerical algorithm used to carry out the characteristic evolution. The properties of the 
close approximation waveforms, are presented in Sec. |V[ 

Notation and Conventions: We use a metric of signature ( — h ++) and a null tetrad 
with normalization /"n^ = —rrf-fha = —1, so that gab = '^{jn(^a'n^b) — ^ant))- We use qab 
to represent the standard unit sphere metric in angular coordinates = ip) and set 
qAB _ q{AqB) ^ where q^^qsc = ^c; with q"^ = (1, i/ siwd). We use q"^ to define the 9 operator 
with the convention 9/ = q'^dAf, for a spin-weight function /. Complex conjugation is 
denoted with a "bar", e.g. 3fJ(/) = (/ + /)/2. The conventions of the present paper result 
in a different form of the perturbation equations from that originally given by Teukolsky. 



II. THE PERTURBATION EQUATIONS 
A. The background Schwarzschild metric in outgoing horizon coordinates 

The Schwarzschild metric in standard coordinates is 

ds'^ = -(1 - + dr^{l - ^)"^ + r^qABdx^dx^. (2.1) 

In outgoing Eddington-Finkelstein coordinates, where £t = t — r* is a null coordinate and 
r* = r + 2Mlog(2^ — 1), the Schwarzschild metric takes the Bondi form 

ds'^ = -(1 - '^)du^ - 2dudr + r'^qABdx'^dx^ . (2.2) 
r 

These coordinates specialize to spherical symmetry the general procedure for constructing a 
Bondi null coordinate system They patch two quadrants of the Kruskal manifold: the 



exterior spacetime quadrant and the quadrant following the initial singularity. 

For the anticipated comparison with a fully nonlinear description of a white hole, it is 
useful to introduce another null coordinate system based upon the affine parameter u = 
-Me""/"^*^ along the ingoing null hypersurface r = 2M that forms the white hole horizon. 
We set M = at the intersection of the black hole and white hole horizons, i.e. at the r = 2M 
bifurcation sphere (corresponding to tt = +oo in Eddington-Finkelstein coordinates). The 
metric then becomes 
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(is^ = -(1 ) — —du^ + dudr + r^qAndx^dx^. (2.3) 

r u 



In addition, we introduce an affine parameter A along the outgoing null geodesies in the r 
direction, with the affine freedom fixed by requiring that A = when r = 2M and that 
g°'^{dau)dh\ = — 1. This implies 

4M(r - 2M) 
A — . 

u 

In these (m. A) coordinates, the metric takes the form 

ds^ = -Wdu^ - 2dud\ + r'^qABdx^dx^, (2.4) 

where 

r = 2M-^ (2.5) 

AM ^ ' 

and 

9X2 

W = — . (2.6) 

These coordinates specialize to spherical symmetry the general procedure for constructing a 
Sachs null coordinate system designed for the double-null initial value problem 0. For this 
reason, they are especially useful for the study of horizons in the nonlinear regime. They 
are also useful in the perturbative regime because they cover the entire Kruskal manifold 
r > with remarkably simple analytic behavior, as first discovered by Israel Since 
g^^ = W = —\^/{2Mr) the hypersurfaces A = const are everywhere spacelike (so that the 
M-direction is spacelike) except on the white hole horizon where A = and the w-direction 
is null. The surfaces u = const are everywhere null. The spacelike surfaces A = const > 
(A = const < 0) can be used as partial Cauchy hypersurfaces to cover the two quadrants 
above (below) A = in the Kruskal manifold. 



B. The Teukolsky equations 

By aligning a complex null tetrad, {l^ , n'^ , m'^ , rh^) with the degenerate principal null 
directions of a Petrov type D background spacetime, Teukolsky was able to express the 
vacuum perturbation equations for the Weyl curvature scalars ipo = Cabcd"'m^l''rn'^ (of spin- 
weight s = 2) and ?/'4 = Cabcd^^'^^n'^^'^ (of spin-weight s = —2) of the Newman-Penrose 
(NP) formalism as the simple wave equations 



[(D + 3e - e + 4p + p) (A + 47 - /i) - (5 - vr + a + 3/3 + 4r) (5 - TT + 4a) - 3^2] i^o = 0, (2.7) 

[(A - 37 + 7 - 4p - p) (D - 4e + p) - (5 + f - /3 - 3a - 47r) (5 + r - 4/3) - 3^2] ^4 = 0. (2.8) 

Here the spin coefficients a = {la-bn°'fh^ — ma-hrh°-m^)/2, (3 = (/a^bn^m^ — ma-hm"'m^)/2, 
7 = iL;bn°-n^ - ma;bm"-n'')/2, e = {h-bn^l^ - nia-bfriH^) /2, r = la-bm^n^, tt = -na-bm°-l^, 
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p = la-^bm^fn}' and p = —na-fiifi""w}' are computed using the background geometry and the 
directional derivatives are D = l°'da, A = n^-da, S = m°'da and 5 = fff'da- The Weyl scalars 
ipQ or '?/'4 are first order quantities in perturbation theory while il)2 = Cabcdl°"<T^^'>TT''^n'^ = —M/r^ 
is a zeroth order curvature quantity. 

This formulation has several advantageous features: (i) It is a completely first order 
gauge invariant description (i.e. the perturbative Weyl scalars tpo or ?/'4 are invariant not only 
under infinitesimal coordinate transformations but also under null rotations of the tetrad); 
(ii) It does not rely on any frequency or multipole decomposition (i.e. the above equations 
can be directly integrated in the time domain); (iii) The Weyl scalars are objects defined in 
the full nonlinear theory, where ipQ can be prescribed as constraint-free data on an outgoing 
null hypersurface and ip4 as constraint-free data on an ingoing null hypersurface. 

Since Eq's. ( |2.7| )-( pl8D are expressed in a covariant form, it is straightforward to write 
them explicitly in any background coordinate system. Specializing thus to the Schwarzschild 
background metric in the Israel coordinates introduced in Eq. (^.4|), we choose a null tetrad 



n 



m 





i-r 




W d 


i-r- 


Y^dx 


1 .a 





[0,1,0,0] 



where 



' di) sin ^ dip ' 



[0,0, g^] . 



(2.9) 



(2.10) 



a 



Correspondingly, the only non-vanishing background NP quantities are, 
cot{^) (r + 2M)A 



2^27 



7 



8Mr2 



D- — A - — — 
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1 



A 

2^' 



u 



P 



AMr'' 



d id 

+ — 



-\/2r V sin ^ dip 



6 



d 



i d 



V '^'^ sin ■{} dtp 



Substitution of the above NP quantities into Eq's. (| 
tions to 



^2 ^2 

(Lo + ^)^o = 0, and (L4 + ^)^4 = 0, 



8]) reduces the Teukolsky equa- 

(2.12) 



where 



1 



AMr 



AMr 



2Mr2 
7 

2r"2 



AMr 2Mr2 ' ' " 2rM' 

X'di + dudx - -^j^^^u - —Xdx 



'X 



2Mr3 



(2.13) 



:2.111 
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and = —95 is the angular momentum squared operator. In order to treat the radiation 
near X"*" it is also advantageous to consider a boosted tetrad n", m", m"), with 1°" = —V^u 
and h'^ satisfying l°'ha = — 1. We accordingly define boosted Weyl scalars ipo = Cabcd"''fn I'^'m'^ 
and '?/'4 = C abcd'tT-'^'ni'^'iV^'iT^'^ ■ This boosted tetrad is adapted to the affine time u at X"*" rather 
than the affine time u at the horizon. 



C. Spin- weight-zero versions of the Teukolsky equations 



It is useful to convert Eq's. ( p.l2|) - ( p.l3|) into spin-weight-zero equations. Considering 
the commutation relation (55 — 55)?7 = 2s?7 for a spin-weight s function 77, and setting 



^0 



S^<l>o so that SS'j/'o = 9^(6 + 99) $0) the Teukolsky equation for -^o becomes, after 



factoring out an overall 9^, 



AMr 



AMr 



-udn 



^ A(3M - 4:r)dx + ^ 



2Mr 



2rM 



(6 + 99) 
2r2 



$0 = 0. (2.14) 



Similarly, setting ip^ = 9^$4 so that QQip^ = 9^(2 + 99) $4 and after factoring out an overall 
9^, the Teukolsky equation for $4 takes the spin- weight-zero form 



1 
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16M2 + AMr (2 + 99) 



AMr" ' """'^ AMr 2r2 " 2Mr^ 2r^ 

These equations can be re-expressed in terms of the Laplacian 

2A(16M2 - uX) 



$4 = 0. (2.15) 



- uxy 



-dx - 2dudx 



(2.16) 



defined by the metric ds'^ = —Wdu^ — 2dudX induced by the Schwarzschild metric on the 
2-dimensional (k. A) subspace. Discussion of the asymptotic behavior of the Weyl tensor is 
most convenient in terms of the variables Fq = r^$o and F4 = r$4. The spin-weight-zero 
Eq's ( 2.14 ) and (|2.15 ) then reduce to 



where 



and 



{D' + To)Fo = 
{D^ + T^)F, = 0, 



(6M + r)A„ 30M (6 + 99) 
o\ 5 1 5 — 



(2.17) 
(2.18) 



T4 



Mr2 
(6M + r)A 



Mr2 



QM (2 + 99) 



do not contain du terms. 

Asymptotic fiatness requires that Fq and F4 have finite limits at X+. This is consistent 
with the asymptotic forms of Eq's ( p.l7| ) and ( |2.18| ), which asymptote to one-dimensional 
wave equations for solutions whose radial derivative falls off uniformly as 0(l/r^). The limit 
of F4 determines the outgoing gravitational radiation waveform and the limit of Fq is related 
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to the retarded quadrupole moment of the system. More precisely, Mmrtp^ = duN, where 
is the standard definition of the Bondi news function. 
In the hmit A oo, 

(6M + r)A 4 

so that the appearance of this term in Tq and T4 can cause inaccuracy in the numerical 
solution of Eq's ( |2.17| ) and ( |2.18| ) at late Bondi times on X"*" (as —u — * 0). In particular, 



the late time behavior of the radiation waveform, can be more accurately computed by an 
evolution algorithm for the Weyl component ip^. Setting ipQ = S^$o and 1IJ4 = 5^ $4, the 
fields Fq = r^$o and F4 = r$4 satisfy 

{D^ + So)Fo = 0, (2.19) 
{D^ + S4) h = 0, (2.20) 



where 



16M(r-3M) 30M (6 + 99) 4(r - 3M) 30M (6 + 99) 

00 — ^ Ox 5 1 7^ — ^ Or 5 h 



and 



16M(r-3M) 6M (2 + 99) _ 4(r - 3M) 6M (2 + 99) 

'J4 — 5 0\ r- H — Or ^ H 5 . 

In these variables, the deviation of the Teukolsky equations ( |2.19| ) and ( p.20| ) from a 1- 
dimensional wave equation is independent of time at a fixed r. These equations are well 
behaved at X"*", i.e. after a compactified coordinate such as x = 1/r is introduced. Because 
F4 = {u/4M)'^F4, where F4 must be regular throughout the Kruskal manifold (since it 
is constructed with a regular basis), it follows that F4 — as the black hole horizon is 
approached. This facilitates an accurate long term evolution of the waveform using Eq. 

(PO). 

Note also, however, that Fq = (4M/m)^Fo so that Fq is singular on the black hole 
horizon and thus a poor choice of variable for long term evolution. The opposite signs of the 
coefficients of dr in 5*0 and 5*4 are responsible for this behavior, as can be seen by ignoring 
the remaining potential terms and freezing the coefficient of 5^ at r = 2M, so that Eq's 
( Prrgp and (^ reduce to 



{2du -dr- ^)drh = 0, (2.21) 
{2du - dr + l^)drF4 = 0, (2.22) 

in terms of retarded Bondi coordinates. In this approximation, both of these equations 
admit purely outgoing waves F{u). However, an ingoing Fq wave has the exponentially 
singular behavior 

Fq = f{u + 2r)e^«~ 

as an initial pulse /(mq + 2r) approaches the black hole horizon as m — > 00. In contrast, an 
ingoing F4 wave decays exponentially on approach to the black hole. 
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D. Time reflection properties of the Teukolsky equation 



The different forms of Eq's. ( p.l4|) and ( p.l5|) , or Eq's. ( p.l9|) and ( p.20|) , make it clear 



that the Teukolsky equation for the Weyl component ipo in the outgoing null direction /° and 
the Teukolsky equation for the Weyl component in the ingoing null direction ra" are not 
related in a way which makes manifest the time reflection symmetry T of the background 
Schwarzschild geometry, defined by T{t, r, 6, 0) = {—t, r, 6, (p) in Schwarzschild coordinates 
or by T{u, r, 6, 0) = {—v, r, 6, (p) in terms of Bondi retarded and advanced times u = t — r* 
and V = t + r*. The time reflection symmetry could be incorporated explicitly by introducing 
null tetrad vectors = al"" and A^" = (l/a)n" satisfying TL" = — A^". However, the explicit 
form of the required boost. 



2M /2(r-2M) 2M / -2Xu 

makes it clear that such a time symmetric formulation would introduce singular behavior at 
both the black and white hole horizons. 

However, this time symmetric tetrad is useful for formulating the time reflection symme- 
try of solutions of the Teukolsky equations using other tetrads. Let ^4 = CabcdN""m^ N'^m'^ = 
'^{u,r,6,(f)) be a perturbative solution for \l/4. Then the time reflection symmetry implies 
that \l/o = CabcdL°"fn^L'^fn'^ = ^{—v, r, 9, 0) is a perturbative solution for \l/o- This correspon- 
dence maps a retarded solution (no incoming radiation) for \&4 into an advanced solution 
(no outgoing radiation) for \l/o. In terms of the and Weyl components, the solution 
ip4 = ip{u, r, 6, 0) corresponds under time reflection to the solution 



III. NULL DATA FOR A SCHWARZSCHILD PERTURBATION 

The Weyl component ipo can be posed as constraint-free gravitational data on an outgoing 
null hypersurface. Similarly, -04 constitutes constraint-free gravitational data on an ingoing 
null hypersurface. These nonlinear results extend to linearized theory but care must be 
exercised in applying them to the Teukolsky equations. 

In the Cauchy formulation, the Teukolsky equation for -04 is normally chosen in order to 
investigate the outgoing radiation introduced by a perturbation. However, the Hamiltonian 
and momentum constraints prevent the free specification of -04 (or ipo) on a Cauchy hypersur- 
face. Consistent Cauchy data for tp^ must be provided indirectly by a 3-metric and extrinsic 
curvature that solve the constraints. In the Cauchy approach to the close approximation. 



this has been provided by (a limit of) Misner's time symmetric wormhole data P7 |. 

In the double-null formulation of the characteristic initial value problem, data are given 
on a pair of intersecting null hypersurfaces, one outgoing and one ingoing. Null data for the 
Teukolsky equation for -04 can be freely posed on the ingoing null hypersurface but data for 
04 on the outgoing null hypersurface has to obtained indirectly. This can be done, as in the 
Cauchy problem, by first considering consistent metric data in double null coordinate, from 
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which the Weyl data for ^4 can be constructed on both hypersurfaces. This is the method we 
use here to generate two examples of double-null data for the Teukolsky equation: Robinson- 
Trautman perturbations and close approximation data. 



A. Robinson- Trautman perturbations 

The Robinson- Trautman space-times describe an algebraically special but distorted 
and radiating black hole. They provide an important testbed for the computation of a 
general perturbative solution by numerical evolution. In the case of outgoing radiation from 



a black hole of mass M, the metric can be put in the Bondi form |29 



(is2 ^ _ '^)du^ - 2Wdudr - 2rW,Adudx^ + r^qAndx^dx^, (3.1) 

where /C = W^[l — L'^(\ogW)], L"^ is the angular momentum operator and yV{u, x^) satisfies 
the nonlinear equation 

UMduihgW) = W^L^IC. (3.2) 

The outgoing Eddington-Finkelstein form of the Schwarzschild metric Eq. (|2.2|) results 
from setting W = 1. More generally, smooth initial data }V{uo, x^) evolve smoothly to form 
a Schwarzschild black hole horizon. The linearized solutions to the Robinson- Trautman 
equation ( |3.2| ) are obtained by setting W = 1 + and dropping nonlinear terms in 0: 

12M(9s0 = L2(2-L2)0. (3.3) 

For a spherical harmonic perturbation = A{u)Y£fn this leads to the exponential decay 

The corresponding Weyl tensor components for the perturbation are ipo = 0, in agreement 
with the role of 1°' as an algebraically degenerate principal null direction, and 

= 2M^. - 2] |-6M + n^+l) r| (_„^M)--^i"'^"""^''--l8^>1„. (3.4) 

in terms of the affine horizon parameter u. The perturbation vanishes on the black hole 
horizon 7i+ and is singular at I~ . This supplies the data on 7i~ and an outgoing null 
hypersurface u = U- for the evolution of ■04 forward in retarded time. 

For the corresponding time reversed solution, ?/'4 = 0. By applying to the Robinson- 
Trautman perturbations a procedure for mapping an outgoing solution of Einstein's equa- 
tions into an ingoing version [^], we find the solutions 

ds^ = -(C - ^)rf5^ + 2Vdvdr + 2rV Advdx"^ + r'^qABdx^dx^ , (3.5) 
rV 

where v is the advanced Bondi time coordinate, C = V^[l — L'^{logV)] and 

12Md^{\ogV) = -V'^L'^C. (3.6) 
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The linearized solutions obtained by setting V = 1 + satisfy 



l2Md^(t) = -L^{2- L^)(t). 



(3.7) 



For a spherical harmonic perturbation = B{v)Y^rn, this leads to the exponential growth 
B = Bqc" +£-2)/i2M_ rjnj^g corresponding perturbative Weyl tensor component is 



bation vanishes on the white hole horizon and is singular at Nevertheless, it can be 
used to check a (forward in retarded time) evolution algorithm, beginning at a retarded 
time U-, by pasting asymptotically fiat initial null data outside some radius to interior 
Robinson- Trautman data. 



The coordinates introduced by Sachs to formulate the double-null characteristic initial 
value problem are especially useful when one of the null hypersurfaces is a white hole 
horizon . Sachs' coordinate system consists of (i) an affine null coordinate u along 
the generators of ?^~, which foliates 7i~ into cross- sect ions and labels the corresponding 
outgoing null hypersurfaces emanating from the foliation; (ii) angular coordinates which 
are constant both along the generators of 7i~ and along the outgoing rays and (iii) an affine 
parameter A along the outgoing rays normalized by V^mVqA = —1, with A = on In 
the resulting = {u, X,x^) coordinates, the metric takes the form 



In addition, it is useful to set qab = ^"^hAB, where det^hAs) = det^qAs), where qab is the 
unit sphere metric. 

The requirement that 7Y~ be null implies that W = on Ti. and the gauge freedom on 
can be fixed by choosing the shift so that du is tangent to the generators, implying that 
= on 7i~, and by choosing the lapse so that u is an affine parameter, implying that 
dxW = on Ti.^ . In addition to these choices, we fix the affine freedom in u by specifying 
it on a slice S~ of 7i~, which is located at an early time approximating the asymptotic 
equilibrium of the white hole at past time infinity /~. The outgoing null hypersurface 
emanating from S~ then approximates past null infinity X~ . The Schwarzschild metric in 
Israel coordinates (|2.4| ) is obtained in the spherically symmetric case when = and 
hAB = Qab- 

The double-null problem for the close limit of a white hole is posed on the ingoing null 
hypersurface 7i~ and the outgoing null hypersurface JT", which extends to X"*". In order to 
pose the double-null Teukolsky problem for -04 in the perturbative regime, we generate the 
data for 04 from metric data for the nonlinear version of the problem. The metric version 
of the null data consists of the values of the spin- weight-two field J = q^q^hAB on 7i~ and 




B. Close limit initial data 



ds"^ = -{W- gABW'^W^)du^ - 2dudX - 2gABW^dudx^ + gABdx'^dx 



(3.9) 
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J . For a perturbation of a Schwarzschild background, r is given by Eq. (p.5|) . On the 
event horizon 7i~, 



where iiT = Vi + JJ, which reduces in the hnear regime to 

^4 ~ \j^uu. (3.10) 

Off the horizon the expression for 7/^4 is more comphcated and involves W and as well 
as J and r. 

The horizon data for a head-on fission of a white hole, can be obtained from a conformal 
model based upon an ingoing null hypersurface emanating from a prolate spheroid embedded 
in a flat space 0]. Let {r,6,(j)) be standard spherical coordinates for the inertial time 
slices t = constant of Minkowski space. In the close limit, the eccentricity of the spheroid 
vanishes and the Minkowski null hypersurface reduces to the light cone from a sphere t = 0, 
f = a. The perturbation of its conformal null geometry is described, to linear order in the 
eccentricity, by 

T/- sin^ 9 /„ . . X 

j{t,e) = — — , (3.11) 

t — a 

where the relation between t and the affine parameter u on the white hole horizon is 



di f\f-l)' / (5-v/l3)-2f xVv^ 



in terms of 

t — a 



P 



(3.13) 



Here p and a are positive parameters which adjust the affine freedom in the position of the 
Minkowski null cone on the white hole horizon. At early times Eq. ( p.l2| ) implies m ~ t but 
as the Minkowski null cone pinches off at t = a the corresponding affine time on the white 
hole horizon asymptotes to m — 00. In terms of the inverted pair-of-pants picture for a white 
hole fission, the pants legs are mapped to m = 00 so that in the close limit the individual 
white holes are mapped to along the white hole horizon in the Kruskal manifold. The 
details are discussed elsewhere in a treatment of fully nonlinear null data for the general two 
black hole problem |]3T| , |32| . 

Close limit data for J{u, 6) on the white hole horizon is determined by integrating Eq. 
( p.l2| ) and substituting into Eq. ( p.ll| ). These equations allow the rescaling u — > Ku, 
t —>■ Kt, p — > Kp and a —>■ Ka which allow us to set p = 1 without any loss of generality. 
Note, that the rescaling u Ku is equivalent to the time translation isometry u u+const. 
In order to eliminate nonessential parameters, we initiate the integration at the bifurcation 
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sphere u = 0. Then, up to scale, the close data are determined by tq = t|„=o < or in 
terms of J by the parameter 

pJ\u=o 1 

V = f, = - — , (3.14) 

UJ I u=oo 

which is independent of the overall scale freedom J XJ that is factored out in the close 
approximation. The parameter 77 is a scale invariant parameter describing the physical 
properties of the close limit. It determines the yield of the white hole fission. In the time 
reversed scenario of a black hole collision, t] would be related to the inelasticity of the 
collision. No similar parameter seems to appear in the Cauchy description of the close 
approximation in terms of time symmetric Misner data P]. 

The close limit data on the horizon for ?/^4 used in the simulations presented in Sec. 
V Q are obtained by integrating Eq. (p . 1 2|) with a 4th order Runga-Kutta scheme and 



carrying out the substitutions into Eq's ( |3.11| ) and ( ^.10[ ). The data on an early outgoing 



null hypersurface are accurately approximated by setting ip^ = since Eq. (|3.11| ) implies 
ipA = 0{u~^). This approximates the condition on the data that there be no ingoing radiation 
at X^. 



IV. NUMERICAL ALGORITHM 

Before giving the details of the numerical algorithm, we should state the goals we want 
to achieve. For many purposes, it would seem sufficient to evolve the waveform until one can 
read off the first few cycles of quasinormal mode oscillation. This is sufficient in practice to 
compare with a nonlinear evolution, to get the astrophysically relevant part of the waveform, 
to compare quasinormal mode results with those in the literature, etc. Instead we define 
as our criterion of quality the ability to evolve stably and accurately well into the domain 
where the waveform is dominated by a power law, which requires at least a lOOOM of Bondi 
time for our typical data. This turns out to be a rather stringent criterion, which rules out 
a number of numerical approaches which we have tried. In all such approaches, our overall 
strategy is to compactify the outgoing null direction and bring X+ into a finite coordinate 
distance while maintaining regularity of the equations. 

We begin the description of our numerical setup with a discussion of the ingoing null 
geodesies, which forms the basis of our approach to the numerical solution of the Teukolsky 
equation. Then we briefiy describe a few of the algorithms which do not work completely 
satisfactorily for the Teukolsky equation, and explain why this is so. We believe that this 
also provides useful experience for nonlinear studies, where, lacking a stationary background 
geometry, the source of problems may be much less obvious. 

Our present results pertain to the Teukolsky equation for ip^, where describes the 
outgoing radiation through its asymptotic 0{l/r) behavior at The incoming radiation 
at X~ is described by ipQ, which has asymptotic 0(l/r^) behavior at X+. As a result, an 
accurate treatment of the Teukolsky equation for ipQ requires different numerical methods, 
which will be described in a forthcoming paper. 
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A. Ingoing null geodesies 



Null geodesies are fundamental to the design of the numerical algorithm since they are 
the characteristics of the Teukolsky equation. In particular, since we handle the angular part 
of the spin-weight-zero Teukolsky equation by a spherical harmonic decomposition in which 
99 = 1), the relevant characteristics are the radial null geodesies. The outgoing null 

geodesies are automatically built into the characteristic evolution scheme, which is based 
upon a retarded time foliation. The remaining issue is how to effectively incorporate the 
behavior of the ingoing null geodesies into the algorithm. 

Consider first the description of the ingoing null geodesies in a compactified version of 
Israel coordinates {u, x), where x = A/(M+A) so that is located at a; = 1. The analytical 
simplicity of these coordinates is not matched by their numerical convenience. The ingoing 
null geodesies satisfy 

dx ,^ ^2 ^ ('1 ~ ^^^^ 



x: = -(i-^: 



du ' ' 2M 8M(1 -x)-ux 
Near X+ where 1— x = 5<< 1, this reduces to 

db _b 
du u 



(4.1) 



Hence 5, and the separation between neighboring geodesies near T"*", decays linearly with u 
and exponentially with u. Thus evolving for a Bondi time of £t = lOOOM is impossible as 
the geodesies would be within e~^^°5o of each other. 

Not only is the x-coordinate numerically unsatisfactory in the way it compactifies X+, 
the ti-coordinate is also inconvenient in the way it covers the exterior Kruskal quadrant in a 
finite retarded time. This prevents the long term numerical resolution of the ringdown (with 
the characteristic time scale of the lowest quasinormal mode) without using an exponentially 
decaying step size Am. This is simple to fix by using Bondi time u as the time step coordinate. 

The problem with the x-coordinate can be "delayed" by introducing a dynamical grid, 
in which the gridpoints move along the ingoing null geodesies, a strategy that has been 



successful in studies of spherical critical collapse ||3^ . This approach drastically decreases the 
discretization error and is sufficient to evolve for about 100 M, and read off the quasinormal 
ringdown frequency and damping time with good accuracy. However, the Ax intervals 
between neighboring geodesies decrease exponentially and the approach breaks down once 
the separation between neighboring gridpoints falls below the error of the geodesic integrator 
(e.g. machine precision) and a "numerical crossing" of the geodesies effectively occurs. 

Note that in practice the computation of the null geodesies is related to the problem of 
inverting the definition of the tortoise coordinate 

r* = r + 2Mlog(— - 1) 

to compute r, which is also required for our production algorithm as discussed in the next 
section. Both problems are handled numerically by solving the above implicit equation 
iteratively using Newton's method (in terms of the appropriate coordinates). 
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B. Numerical algorithm for outgoing radiation 



The preceding considerations lead us to the following choice of algorithm for a production 
level unigrid code with optimal performance. It is based on a {u, p) coordinate system, where 
p is a radial coordinate, which compactifies X"*", defined implicitly by 

r* = Pq tanp, (4.2) 

with Po an adjustable parameter and — 7r/2 < p < tt/2. The p coordinate allows good 
resolution at all times near both and the white hole horizon, . 

In this coordinate system, we evolve the ith spherical harmonic component of F4 by 
expressing the second order differential equation ( p^.20D as the two coupled equations 

dpF^ = G 

cosV^ ^ / rr-3M) sinpcosp\^ po {r - 2M){{f + i - 2)r + 6M) ~ 

Ou'^ = — Op'^ + ^ 7, Cj — 



2po ^ \ Po J cos^ p 2r^ 

The background mass M can be scaled out of the above equations by the rescaling po Mpo, 
u — s> Mu and r — ^ Mr. In this way, simulations can be carried out with M = 1 without loss 



of generality. Note that rescaling u is independent of the rescaling u Ku (see Sec. pi B|) 
which generates the translation isometry of u. 

The relatively sensitive features of Eq. ( [4.4| ) on a hypersurface of constant u are located 
in the region near r = 2M. This is the main reason why the p coordinate is so useful since 
it concentrates grid points in that region while maintaining a uniform grid spacing. We 
make the choice po = 40M, which gives good resolution throughout the evolution, well past 
the ringdown phase and into the final power law tail. Note that this would not be possible 
with the simpler approach of writing the Teukolsky equation in (a compactified version of) 
double null coordinates (m, w), which gives excellent results until one reaches the power law 
tail. At this stage the dynamics is essentially dominated by the "Schwarzschild potential" , 
which is not well resolved in double null coordinates. 

Equation ( |4.3| ) is solved using a second order accurate integration in p. Equation ( |4.4| ) is 
solved using a second order (in time) Runga-Kutta scheme, with the dpG term evaluated by 
means of second order accurate forward differencing in the interior of the grid, and second 
order accurate central differencing for the point neighboring X+. (The dpG term drops out 
on both IHr and where du is a characteristic direction). 

For a typical choice of initial data the power law tail only sets in after the quasinormal 
oscillations have decayed by more than 10 orders of magnitude. In order for the final tail 
not to be lost in machine error it is necessary to evolve the quasinormal phase in quadruple 
precision. 



V. WAVEFORMS 

In this section, we present computed waveforms for three types of quadrupole data, with 
the background mass scaled to M = 1. The first case, an analytically known Robinson- 
Trautman perturbation, is used to establish second order convergence of the numerical algo- 
rithm. The Robinson- Trautman waveform decays as a pure exponential. The second case, 
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a pulse of compact support emanating from the white hole horizon, serves to monitor the 
ability of the code to track many cycles of the quasinormal ringdown of a generic radiation 
tail. The final case is the close limit waveform from a white hole fission. 

A. Robinson- Trautman testbed 

An ^ = 2 Robinson- Trautman perturbation is determined by the spin-weight-0 field 

F4 = (5.1) 

r 

We use this to perform a convergence test of the code by evolving from -u = to tt = 5 and 
then examining the £00 norm of the error E = \\F4^,numeric ~ F^^analyticW versus grid 
spacing at £t = 5, at which time the signal has decreased by a factor of e~^^. The convergence 
plot of the error given in Fig. determines a slope of 2.0041, in excellent agreement with 
the theoretical second order accuracy of the algorithm. 

B. Ringdown from a compact pulse 

We simulate the evolution of an £ = 2 quadrupole pulse of compact support emerging 
from the white hole horizon Ti.^ . The pulse consists of a single peak of the form 

F4{U,X = 0) = {{U - Umin){Umax , Umin < U < Umax, (5.2) 

with F4{u,x = 0) = outside this interval. Figure |^ shows the waveform on X"^ obtained 
by evolving with initial data F4 = on an outgoing hypersurface preceding the pulse. In 
the simulation used to produce the waveform at X"^, we choose Umin = —50 and Umax = 
and evolve from u = — 60M to u = 2000M . The simulation was performed in two steps. 
The first step, in the interval from u = — 60M until u = 250M, was performed in quadruple 
precision in order not to lose the final tail in roundoff error. The second step, in the interval 
from u = 250M until u = 2000M, was performed in double precision. 

Figure ^ is a logarithmic plot of the absolute value of the waveform versus u for the same 
data. It covers the period from the onset of quasinormal ringdown to the onset of the final 
tail decay. The logarithmic plot clearly demonstrates the exponential decay and shows a fit 
to a quasinormal decay. The lowest quasinormal mode for a gravitational perturbation of the 
Schwarzschild metric has the theoretical form / ~ sin(.373672'u) exp(— .0889625n) |Q. The 
fit of the computed waveform to a quasinormal decay is F4 ~ sin(.373668'u) exp(— .088951u), 
in excellent agreement with the expected theoretical form. The corresponding fit of 
the close approximation waveform given in Sec. [VQ yields the quasinormal dependence 
F4 ~ sin(. 3736735-11) exp(— .0889575-0). A conservative comparison of these two calculations 
indicates a quasinormal dependence F4 ~ sin(. 37367m) exp(—.08895tt), with the numerical 
uncertainty in the last digit. 

Figure ^ shows a log-log plot of the late time behavior of the waveform and the final tail. 
The measured slope of the tail indicates a power law decay , with the power varying from 
F4 (X -u~^'^^ near the beginning of the tail to F4 oc u'^-^^ near the end. 
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C. The close approximation waveform 



As discussed in Sec. [ill the effective parameter space for the head-on close approxi- 
mation data can be reduced to the single scale invariant parameter rj controlling the fission 
yield. In the simulations presented here we set the scale dependent parameter p = 1. Thus, 



in accord with the discussion in Sec. IIIB, we identify t — a = f and prescribe the horizon 



data used in the simulations in the normalized form 

F4 = -(A9,)(A9,)4, (5.3) 

T 

after factoring out the i = 2 angular dependence. 

The time dependence of the close approximation data is quite mild when expressed as a 
function of f, as in Eq. ( |5.3|) . However, the relationship ( p. 121) can cause the dependence on 



u to be quite sharp. There is a transition region where the behavior of f(u) changes from 
the asymptotic form df/du ^ 1 as f — —oo to df/du — > as f ^ 0. For large values of 
the parameter 77, this produces sharply pulse shaped data, as described below. 

Figure |^ plots f versus u for 77 = 7060, 1410, 364, and 84.3. The plots reveal a relatively 
sharp transition in the slope. This transition region is translated in the negative w-direction 
as 7] increases. For sufficiently small 77 the transition occurs at m > 0, in the region of the 
white hole horizon which does not affect the exterior spacetime. 

The location of the transition region affects the nature of the horizon data. Figure ^ 
shows F/^{u)\'yi- for the above values of rj. The value of 77 only changes the position of the 
transition region, not its width. Hence a change in 77 translates the horizon data F^{u)\-}-i- 
but does not change its shape. 

The horizon data for Fi{u) has a more complicated dependence on rj due to the expo- 
nential relationship between u and u and the extra factor of introduced by the change in 
tetrad. This is of physical importance since it is F4^{u) which is the observed waveform at 
X"*". The factor of v? suppresses pulses centered at smaller \u\ compared to those centered 
at larger \u\ and forces the resulting pulse to vanish at m = 0. The relation between u and 
u varies from an exponentially increasing blueshift at large negative u to an exponentially 
increasing redshift at m = 0. This has the effect of compressing pulses centered at more 
negative u compared to those centered at less negative u. These effects combine to pro- 
duce successively broader pulses for successively smaller rj. However, once rj is sufficiently 
small, the transition region is located at m > and rj does not affect the shape of F4(£t)|-^-, 
although it affects its overall amplitude. In that case f u + tq for —00 < u < 0, and 
Fi\H~- oc l/(M + fo)'^, where tq = — 1/?7. As a result, modulo a constant overall multiplicative 
factor and a constant shift in £t. 

Figures ^ and || show F^{u)\-yi- and F4('u)|-h-, respectively, for rj = 7060, 1410, 364, and 
84.3. Figure |] shows Fi\n- versus u for small rj, with the amplitude and position of the 
peaks adjusted so that they overlap. Except for the overall amplitude, there is no significant 
effect on the data even when rj is reduced by 3 orders of magnitude. For small values of 
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T], Fi(u)\-^~ scales linearly with r] in accord with Eq. ( ^.4| ), whereas for large r] it scales 
quadratically, as evident from the renormalizations in Fig's ^ and |^. 

We test the convergence of the waveform at X"*" by evolving this close approximation 
data with increasingly larger grids containing 1001, 2001, and 4001 points. We define 6yi 
to be the difference between the waveforms obtained using 1001 and 2001 points, and Sy2 
the difference between using 2001 and 4001 points. (We consider only points common to all 
three grids.) Second order convergence requires that 6yi = 4(5?/2- For these grid sizes. Figure 



rO| shows that 5yi and 4(5?/2 overlap confirming that the code is second order convergent 



throughout the quasinormal ringdown phase. 



Figure |TT| shows a series of waveforms produced on X"*" obtained by evolving the close 
approximation data for rj = 7060, 1410, 364 and 4.39. The waveforms have been translated 
with respect to each other and normalized to unit amplitude for purpose of comparison. The 
plots show the waveforms from the initial time up to (roughly) the onset of quasinormal 
decay. 

Figure [l^ shows a log plot of the waveform produced for rj = 158. The fit of the 
exponentially damped section is 

F4 oc e-™^^" sin(.3736735{i), (5.5) 

which matches the the theoretical form for the lowest quasinormal mode to five digits (in 
the frequency). 

Figure |13| shows the late time tail of the waveform. The measured slope of the tail 
indicates a power law decay of the approximate form F4 oc u~^'^ near the beginning of the 
tail to F4 oc u~^'^ near the end, very similar to the behavior of the tail for the compact pulse 
described in Sec. |V lj| . These results suggest a final integer power law tail F4 oc u~^. For an 
£ = 2 quadrupole wave, this is the same ■u~(2^+2) integer power law originally predicted by 
Price for the decay of an initially static multipole. A rigorous mathematical treatment 



of power law tails has not yet been given |^ and it would be particularly interesting to 



reexamine the theory in the context of our boundary conditions. 



VI. DISCUSSION 

Our results establish the capability of characteristic evolution of the Teukolsky equation 
to determine an accurate advanced solution for the head-on collision of black holes in the 
close approximation. In subsequent work, we will extend these results to determine the 
physically more appropriate retarded solution. In the fully nonlinear regime, the conformal 
horizon model for suppling binary black hole data, combined with an existing characteristic 
evolution code, offers a new way to calculate the merger-ringdown waveform from coalescing 
black holes. Because this is an unexplored area of binary black hole physics, these perturba- 
tive studies of the head-on collision will provide a preliminary physical check on extending 
the work to the nonlinear and nonaxisymmetric case, where inspiraling black holes can be 
treated. 
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FIG. 2. Waveform F4{u) at X+ produced by a single pulse emerging from the white hole horizon. 
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FIG. 6. Close approximation data F4(u) on 7i for 4 values of ij. 
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FIG. 7. Close approximation data F^iu) on n~ for 77 = 7060, 1410, 364 and 84.3, with the 
amphtudes renormahzed by the relative factors of 1, 24.03, 305.9 and 3223, respectively. 
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FIG. 8. Close approximation data F4{u) on H~ for 77 = 7060, 1410, 364 and 84.3, with the 
amphtudes renormahzed by the relative factors of 1, 24.03, 305.9 and 3223, respectively. 
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FIG. 11. Close approximation waveforms F4{u) on X+ for 77 = 7060, 1400, 368 and 84.3, with 
the ampHtudes renormahzed by the relative factors of 1, 37.9, 544 and 11200, respectively. 
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FIG. 13. Close approximation waveform F^^u): Late time power law tail for rj = 158. 
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